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Abstract 


This report presents a unified method for subsonic and supersonic jet analysis 
using the three-dimensional Navier-Stokes code PAB3D. The Navier-Stokes code was 
used to obtain solutions for axisymmetric jets with on-design operating conditions at 
Mach numbers ranging from 0.6 to 3.0 , supersonic jets containing weak shocks and 
Mach disks , and supersonic jets with nonaxis ymmetric nozzle exit geometries . This 
report discusses computational methods, code implementation, computed results , and 
comparisons with available experimental data. Very good agreement is shown 
between the numerical solutions and available experimental data over a wide range 
of operating conditions. The Navier-Stokes method using the standard J ones-Launder 
two-equation k-E turbulence model can accurately predict jet flow , and such predic- 
tions are made without any modification to the published constants for the turbulence 
model. 


Introduction 

Knowledge of jet mixing aerodynamics is vital to 
several areas of commercial and military aircraft design, 
such as jet propulsion efficiency, propulsion integration, 
aeroacoustics, and jet interference with aircraft structure. 
Initial jet flow conditions are determined by nozzle exit 
pressure, temperature, Mach number, and nozzle geome- 
try. Once the flow leaves the jet nozzle, the jet flow 
becomes a free shear layer. The action of turbulence 
dominates flow developments farther downstream. As 
such, jet flow properties are difficult to measure or pre- 
dict analytically. 

When Prandtl introduced his mixing length hypo- 
thesis for turbulent flows, a brief analysis of a fully 
mixed jet was given as an example. Early analyses of jet 
mixing behavior were based mainly on this mixing 
length hypothesis and one-dimensional momentum the- 
ory. (See refs. 1-4.) Mean flow properties derived from 
these analytical models compare well with experimental 
measurements of jets at low subsonic speeds. However, 
data from jet flow measurements in the high subsonic 
and supersonic speed ranges (ref. 5) indicate signifi- 
cant departure from the results obtained by using one- 
dimensional momentum theory. 

Jet flow contains a rich combination of flow inter- 
actions and flow physics. These combinations include 
turbulent mixing and compressibility effects such as 
isentropic expansion and shock. Other factors may 
include chemical reactions or shear layer instability. 

Subsonic jet flow features are relatively simple. The 
main variable in the flow is shear layer development 
along the stream wise direction. The static pressure value 
is almost constant with the ambient pressure. In the 
absence of a pressure gradient, no significant in viscid 
flow feature will appear in a subsonic jet. According to 
reported experimental measurements, all turbulent axi- 
symmetric subsonic jets below Mach 0.6 are similar if 


the flow variables are normalized by jet density and noz- 
zle exit velocity. 

On the other hand, supersonic jet flow features can 
be very complex. Because of the supersonic nozzle exit 
Mach number, jet exit pressure can differ from ambient 
pressure. This pressure difference between the jet and the 
ambient fluid must be resolved locally either across an 
oblique shock, by a prominent streamline curvature at the 
jet boundary, or by a Mach disk inside the jet. In addi- 
tion, shocks formed near the nozzle exit may reflect 
repeatedly at the sonic line in the shear layer. Although 
the convected turbulence interacts with shocks in the jet, 
the position of the reflected shock depends on the loca- 
tion of the sonic line in the turbulent shear layer. Such 
interdependence of flow interactions can become very 
complex. 

Earlier jet flow analysis codes, with or without 
chemical interactions included, were formulated with 
simplified assumptions of the Navier-Stokes equations 
and the turbulence model to provide the best jet flow 
simulation within modest limits of computing resources 
available during this time. Analytical methods and simu- 
lation codes developed by this approach have been suc- 
cessfully applied to problems in air-breathing engine 
development, acoustics, and rocket propulsion. (See 
refs. 6-12.) However, there are some drawbacks to this 
approach. First, simplified assumptions are often difficult 
to justify. Second, application of the simplified formula- 
tions is limited to jet flow simulation. The formulations 
are difficult to integrate with computational codes for air- 
frame aerodynamics when performing propulsion air- 
frame integration analysis. It is preferable in such cases 
to perform the analysis with the three-dimensional 
Navier-Stokes equations without empirical assumptions 
for jet flow alone. 

For general use of jet flow simulation, some basic 
requirements must be met. The Navier-Stokes code 
should be upwind biased to capture internal shocks and 



other jet flow discontinuities. The code should also be 
fully three-dimensional in space because the relations 
between turbulent kinetic energy and Reynolds stresses 
are basically three-dimensional. The turbulence model 
should be capable of providing a time scale and a 
consistent description of the production and transport 
properties of turbulent kinetic energy. Therefore, a two- 
equation turbulence closure model is required. 

Many upwind-biased three-dimensional Navier- 
Stokes codes are available that meet the jet flow simula- 
tion requirement. However, the availability of codes with 
a robust two-equation turbulence mode in this class is 
limited. In this report, the PAB3D code is used for all jet 
flow computations. The purpose of this report is to show 
the feasiblity of establishing a unified method for sub- 
sonic and supersonic jet analysis with a general purpose 
three-dimensional Navier-Stokes code. 

The PAB3D code is developed to obtain numerical 
solutions to the Reynolds averaged Navier-Stokes equa- 
tions in three-dimensional spatial domain. The main 
solver algorithm is the upwind Roe scheme, for which 
the numerical dissipation is small. The Jones-Launder 
(ref. 13) two-equation k-e turbulence closure model is 
used to compute the turbulent stresses in the flow. This 
approach is chosen for jet flow analysis because it is con- 
sistent in tracking production and transport properties of 
turbulence kinetic energy and dissipation scale length in 
the shear flow. In the Jones-Launder k-E turbulence 
model, several empirical constants are required. Only the 
published values for these constants are implemented in 
the PAB3D code. These values are fixed for all computa- 
tional applications of the PAB3D code. 

This report describes the mathematical formulation 
of governing equations, the turbulence model, and the 
adaptive grid generation algorithm, along with the 
numerical implementation of each. The adaptive grid 
generation algorithm is designed especially for nonaxi- 
symmetric jet computations. 

Several categories of jet flow computations are 
described separately in the section “Results and Discus- 
sion.” The first category describes axisymmetric jets 
operating at on-design exit conditions so that the jet exit 
pressure matches the ambient static pressure. Results are 
obtained for jet exit Mach numbers ranging from 0.6 
to 3.0. Computed velocity and turbulence intensity distri- 
butions in the jet are compared with experimental data. 
The second category presents results for supersonic jets 
with internal weak shocks. The discussion includes com- 
puted results for jet exit pressures above and below the 
ambient static pressure to show characteristics of the 
shock-containing supersonic jets. Computed results are 
compared with available experimental data. The third 
category of computed cases is axisymmetric supersonic 


jets with embedded Mach disks. Flow conditions for 
these jets are the result of a supersonic jet nozzle operat- 
ing at pressures far from nozzle design value. One partic- 
ular case details a Mach 1.5 nozzle operating at a nozzle 
pressure ratio 3.15 times greater than the design value for 
this nozzle. The last category includes supersonic jets 
with nonaxi symmetric initial cross sections. Shear layer 
development of these jets is very different from a typical 
axisymmetric jet because of added geometrical degrees 
of freedom. The development of elliptic, rectangular, and 
square jets operating at the same exit pressure and Mach 
number is compared. 


Symbols 


a local speed of sound 

e internal energy per unit mass 

Cj, C 2 , constants in two-equation turbulence 
model 


F, G, H inviscid flux components in Navier-Stokes 
equations 

F v , G v , H v viscous flux components in Navier-Stokes 
equations 

A A A 

F, G, H total flux vectors (inviscid plus viscous) in 
Navier-Stokes equations 

fbfl'fn monitoring function for grid adaptation 

Uh k grid index in J;-, r|-, ^-directions 

J Jacobian of coordinate transformation 

k turbulent kinetic energy 

L c jet potential core length 

^1 shock cell length measured from nozzle 

exit to first shock intersection at jet 
centerline 

M Mach number 


NPR 

n 

P 

P 

Pe 

Po 

Pt 

Q 

R 


nozzle pressure ratio, — 

Po 

distance in a direction normal to a solid 
wall 

production term for turbulent kinetic 
energy 

pressure 

jet exit static pressure 
ambient static pressure 
jet total pressure 

conservative variable vector in Navier- 
Stokes equations 

jet exit radius or area equivalent radius 
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S source function in Navier-Stokes 

equations 

T temperature 

t time 

\J e jet exit velocity 

U c jet centerline velocity 

m, v, w velocity components in x- y y- t and 

z-directions 

u'j components of turbulence velocity 

fluctuation 

u fms root-mean-square value of turbulence 

velocity fluctuation, Ju '■ u'- 
u t turbulence velocity fluctuation 

W shock cell length measured from nozzle lip 

to position of first shock reflection in shear 
layer 

x , y, z spatial coordinates 

r compressibility correction factor 

5 ij Kronecker delta 

e turbulent kinetic energy dissipation rate 

p dynamic coefficient of viscosity 

v kinematic coefficient of viscosity 

rj, £ generalized coordinate as function of x, y, 

and z 


turbulence closure models for practical applications. 
Because one of the dominant factors governing jet 
dynamics is turbulent shear layer mixing, the turbulence 
closure model is essential for realistic jet flow simulation 
when using Navier-Stokes methods. The Jones-Launder 
(ref 13) two-equation k-Z turbulence model is used in 
this study. The Navier-Stokes equations and the mathe- 
matical representation of the two-equation turbulence 
model are described briefly in separate subsections of 
this report. 

For computation of nonaxisymmetric jet flows, a 
special requirement in grid generation arises. High grid 
density is required for regions occupied by the shear 
layer and the embedded shock so that high gradients of 
mean flow and turbulence quantities can be accurately 
represented in the numerical solution. However, the posi- 
tion of the shear layer and the shock positions of a non- 
axisymmetric jet are not known in advance. This special 
requirement can be met by using an adaptive grid. The 
analytical basis for an adaptive grid is described in the 
report section “Grid Adaptation Strategy” following 
discussions of Navier-Stokes equations and the Jones- 
Launder k-z turbulence model. 

Navier-Stokes Equations 

The mass, momentum, and energy conservation 
equations of the Reynolds averaged Navier-Stokes equa- 
tions can be written in terms of generalized coordinates 
and in a conservative form as follows: 


p density 

o e , G* constants in two-equation turbulence 

model 

T shear and normal stress components 

Subscripts and superscripts: 
e jet exit condition 

k turbulent kinetic energy 

L laminar quantities 

o free-stream condition 

T turbulence related quantities 

v viscous component of flux vectors 

£ turbulent energy dissipation 

Governing Equations 

The governing equations of the Reynolds averaged 
Navier-Stokes formulation include the conservation 
equations for mass, momentum, and energy and the 
equation of state for gas. In this study, the perfect gas law 
is chosen to represent the properties of air. For a turbu- 
lent flow, the Reynolds stresses can be represented by 


22 22 + 2 = c 
dt ^ d; 3r| 


( 1 ) 


where t, rj, and £ are the independent variable for time 
and the general curvilinear coordinates in the grid 
domain, Q is the conservative flow variable vector (p, 
p«, pv, pw, pe) in generalized coordinates, F, G, H are 
the total generalized flux vectors including inviscid and 
viscous components, and the source term S is zero for the 
Navier-Stokes equations in this form. This equation is 
introduced here mainly to indicate the relationship 
between the basic Navier-Stokes equations and the two- 
equation turbulence model equations. Reference 14 pre- 
sents details of the Navier-Stokes equations as applied in 
the PAB3D code. A simplified form of the Navier-Stokes 
equations which omits all the streamwise derivatives of 
the Reynolds stresses is used in the PAB3D code. Omis- 
sion of these terms is done for computational economy 
and does not introduce significant computation error. 
The remaining cross stream derivatives are numerically 
implemented at several levels in PAB3D. The thin layer 
Navier-Stokes approximation is one option for the user. 
This study uses the option of uncoupled Reynolds stress 
derivatives in two directions. 
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Jones- Launder Two-Equation Turbulence Model 

The Jones-Launder formulation for the two-equation 
turbulence model uses the turbulent kinetic energy k and 
the dissipation rate e as the principal variables. This 
study uses the expanded three-dimensional form (ref. 15) 
of the original Jones-Launder model. This modified for- 
mulation is fully three-dimensional, and the governing 
equations are written in a conservative form as general- 
ized coordinates. The governing equations can be cast in 
the same form as the Navier-Stokes equations, where 


pe 

p* 


f = \ P«e 

p uk 


pve 

pv/: 


H = \ P w£ 

pH’* 


de 1 


3e 



** E d;r 

G v = - 

N* 

H„ = - 

^ z dl 

dk 


dk 

^ Ty 

V 

dk 


3w 


dv 


p = x Lt + T --- + T - 


'dx yydy zz dz 


l T 

xy 


^ du + dv'' 
y dv dx. 


+ x. 


yz 


dv dw 

. T T 

dw du 

dz + dy J 

+ T zx 

^3* dz } 


S = 


S e = C l P l~C 2 pf 


e-2v 


djk 

v" 5 "/ 




+ 


s k = p-p{ \ +r)e + L k 

Here, P is the full three-dimensional production term 
defined as 


P = + , T r^vv T 

xx dx yy dy zz dz xy 


t'du dv ^ 

y dy + dx } 


+ T 


yz 


^dv dw^ 
dz dy 


+ T i 


( dw du ^ 
3.x dz j 


or is expanded to 

o ..rlYdw . dv y 



dw du¥ 


3* 3z 


+ 2 ^ ^ 


r dw ^ 
K dz J 


du dv dw 
dx dy dz 



du dv dw 
dx + dy + dz 


where 

V T = C^Pj p e = + 

b c e 


p k = 

V L 



c ,= 

0.09 

1.44 

C 2 = 1.92 

Q 

II 


Y \ 


T T. = uT 

du i duj 


U r 

dx . 3x. 

\ J 





In the definitions of S E and S h the terms L e and L k are 
near- wall effects which are not important to free jet cal- 
culations, and d/dn denotes derivatives in a direction 
normal to the solid wall boundary. However, these terms 
are included in the PAB3D code. The function Y is the 
compressibility correction function. Several corrections 
have been developed by different authors. Among the 
widely used compressibility correction functions T are 
those proposed by Sarkar et al. (ref. 16) and by Wilcox 
(ref. 17). 

Sarkar model (ref. 16): 

r = ( 2 ) 

Wilcox model (ref. 17): 

r = (Af2 - M 2 0 )H{M T - m t q ) (3) 

where H is the Heaviside function, M T = Jk/p/a is 
the local turbulence Mach number, a is the local speed of 
sound, and Mj 0 is a cutoff turbulence Mach number. 
The commonly accepted value M r 0 - 0.25 is used in 
the PAB3D code. The compressibility correction factor 
is required when the local flow Mach number is greater 
than 1.0. In the Sarkar model, the compressibility correc- 
tion is activated everywhere in the flow field when 
applied for a given computation. The Wilcox model is a 
modification of the Sarkar model so that Y is nonzero 
only for values of M T greater than M J o . This condi- 
tion implies that compressibility correction is activated 
for local flow Mach numbers near or greater than 1 , with 
no correction otherwise. 

Grid Adaption Strategy 

For an accurate representation of the flow field, suf- 
ficient grid density must be provided in the mixing 
region. Unlike an axisymmetric jet, the nonaxisymmetric 
jet is not self-similar and can evolve in dramatically dif- 
ferent fashion in different sectors of the jet cross section. 
Because the position of the shear layer is not known in 
advance, a large number of predetermined grid points 
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are required to provide high density coverage of the 
three-dimensional space if a fixed grid is used for the 
computations. An alternative is to provide high grid den- 
sity in the appropriate locations by using an adaptive grid 
strategy. 

For jet plume analysis, high grid density is required 
in high velocity gradient regions in the shear layer and 
in high pressure gradient regions near shock fronts. The 
number of grid points in each direction of a structured 
grid is fixed. Local grid density can be varied by redis- 
tributing the available grid points in the computational 
domain to match selected flow characteristics such as 
pressure and velocity gradients. Various methods can be 
used to redistribute grid density according to given multi- 
ple functional requirements. In this study, the equi- 
distribution principle and the alternate direction grid 
adaption method published by Eiseman et al. (refs. 15, 
18, and 19) are chosen as the basis for adaptive grid 
implementation in the PAB3D code. 

In the equidistribution approach, a monitoring func- 
tion which governs grid density over the computation 
domain is defined. The monitoring function can be geo- 
metrically represented as a hypersurface in a space with 
dimensions that are one higher than the spatial dimen- 
sions of the computational domain. The process of grid 
adaption begins by constructing a uniform mesh over 
the monitoring surface. For a one-dimensional case, the 
monitoring surface is a curve over the linear spatial 
domain. Equidistribution is simply a uniform distribution 
of points at equal arc distances on the entire length of the 
monitoring curve. When this distribution of points is pro- 
jected back to the one-dimensional baseline in the physi- 
cal domain as adapted grid points, the grid density is 
proportional to the gradient of the monitoring function. 
For a two-dimensional grid domain, the monitoring sur- 
face is a curved surface in three dimensions over the two- 
dimensional physical space. The equidistribution process 
involves constructing a mesh system over the curved 
monitoring surface so that all the grid cells enclose 
approximately equal surface areas. Once the equidistri- 
bution construction is complete, the mesh pattern on the 
monitoring surface is projected onto the original physical 
domain. Similar to the one-dimensional case, high grid 
densities are again obtained in regions where the moni- 
toring function has high gradients. 

If the monitoring surface is assumed to represent 
mountains and valleys in a landscape, the previously 
mentioned process is similar to making a contour map of 
this landscape. Steep slopes in the landscape are natu- 
rally represented by tightly packed contour lines on the 
map, which is a horizontal projection of the original 
three-dimensional surface. Visualization of the monitor- 
ing surface can be difficult for a three-dimensional spa- 


tial domain. However, the algebra and the geometry for 
the adaption process remain unchanged. In addition to 
equidistribution of arc lengths or areas on the hyper- 
surface, normal curvature of the monitoring surface can 
also be used as a weighting function to provide additional 
mesh density control. 

The alternative direction adaption proposed by Eise- 
man simplifies the equidistribution procedure by per- 
forming arc length equidistribution on the monitoring 
surface along each family of coordinate lines. If cell 
skewness remains within reasonable limits, equidistribu- 
tion of all sides of the grid cells will also distribute the 
cell area or volume into approximately equal sizes. How- 
ever, orthogonality is not enforced in this procedure. The 
degree of grid concentration for given values of the gra- 
dients of a monitoring function is controlled by a propor- 
tional constant. Since orthogonality is not enforced in the 
alternate direction equidistribution procedure, excessive 
cell skewness and cell collapse can occur if the propor- 
tional constant is given too high a value. 

For grid adaption to more than one flow quantity, 
multiple monitoring functions can be used. A simple 
approach is to combine all monitoring functions as a 
single weighted sum. The approach of Eisman and 
Brockelie (ref. 20) treats each monitoring function as an 
additional geometrical dimension (which is orthogonal to 
all previous dimensions). In this approach, grid features 
represented in each monitoring function remain distinct. 
The differential arc length element can be given as 

ds x = 1 + w(s 0 )J 1 + |grad(/j)| 2 + |grad(/ 2 )| 2 + ~ ds Q (4) 

where ds 0 is the arc length in the physical or grid domain, 
ds x is the arc length on the monitoring surface, grad 
denotes a component of the gradient in the tangential 
direction of the coordinate curves, and w(j 0 ) is an 
optional weighting function which is proportional to the 
curvature of the monitoring surface. 

A modified approach called the sequential adaption 
method is used in this paper. Assuming there are 
N monitoring functions, the monitoring functions are 
applied sequentially. After each step of adaption, the 
mesh on the previous monitoring surface is treated as a 
“stretched” uniform mesh to support the next monitoring 
function. The arc length increments on each of the moni- 
toring surfaces can be written as 

ds„ = 1 + **'„(•'„ _ i)7 ! + | s rad (/ n )| 2 ds n- l ("= 1-2, 

(5) 

Once the adaption process is completed over the last 
monitoring function, the mesh coordinates are projected 
sequentially back to all previous base surfaces. The 
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last one is the physical space where an adapted grid is 
established. 

Where the curvature weighting functions wj, w n 
are zero, the sequential adaption method and the vector 
monitoring function method are mathematically the 
same. Only their geometrical interpretation and computa- 
tional implementation are different. 

Computational Methods 

The simplified Reynolds averaged Navier-Stokes 
equations and the Jones-Launder k-e turbulence model 
are implemented in the PAB3D code for general fluid 
dynamics analysis in three-dimensional space. Distinc- 
tive features of this code include provisions to accept a 
multiblock grid with patched interface, compact memory 
requirement, and solver options. In particular, a space- 
marching solver with adaptive grid capability is available 
for jet flow computation when flow conditions meet the 
space-marching scheme criterion. For such cases, the 
space-marching solution accuracy is indistinguishable 
from accuracy obtained by using the time-marching 
solver algorithm. The space-marching procedure can 
complete a converged solution in approximately one- 
twentieth the computer time required by the time- 
marching solver to obtain a solution for the same flow 
conditions. 

Solver Algorithm 

Three numerical schemes have been implemented in 
the PAB3D code as solvers for the Navier-Stokes equa- 
tions: the Van Leer flux-vector splitting scheme, the Roe 
flux-difference splitting scheme, and a space-marching 
scheme that is a modified version of the basic Roe 
scheme. These schemes are implicit, upwind, and con- 
structed by using the finite volume approach. Only the 

inviscid portion of the flux vectors F, G, and H is sub- 
jected to the splitting and upwind procedures. The diffu- 
sion terms of the Navier-Stokes equations are centrally 
differenced. Reference 14 details mathematical descrip- 
tion of these schemes. 

The flux-vector scheme and the flux-difference split- 
ting scheme are used in all three computational direc- 
tions. An updated solution at each iteration is obtained by 
using an implicit procedure in the mesh T], ^-planes at 
constant values of ^ and a relaxation procedure in the 
^-direction consisting of a forward and a backward 
sweep. This particular implementation strategy has an 
advantage for computational efficiency. Since the met- 
rics for the implicit procedure are required for only up to 
three planes, the metric constants are recomputed one 
plane at a time at the advancing front of the prevalent 
sweep direction instead of being stored for the entire grid 


domain. Moderate or large mesh sizes require an average 
of only 22 words of memory for each gird point. This 
highly efficient use of computer memory is obtained at a 
modest cost of approximately a 3-percent increase in 
computer time per iteration. The overall computer time 
requirement per iteration per grid point is similar to other 
codes of this type. 

For time-marching solutions using the Van Leer or 
the Roe scheme, each iteration count consists of either a 
forward or a backward sweep in the ^-direction with one 
step of implicit update of the solution in each of the cross 
planes. The inviscid terms in the Navier-Stokes equa- 
tions in the Roe scheme are cast in the form of an 
approximate Riemann problem. The interface flux in the 
streamwise direction is determined by separate terms, 
depending on the quantities on the left (upstream) and the 
right (downstream) sides of the interface. For a fully 
supersonic flow, the information can travel only in the 
flow direction. Such information is carried by the terms 
representing upstream dependence. The terms which 
carry the downstream dependence can be ignored with- 
out introducing significant flow solution error. This state 
of information transfer in the Roe scheme solver is true 
for a broad category of subsonic and supersonic jet flows 
where the streamwise pressure gradient is small. By 
ignoring the downstream dependence terms in the Roe 
scheme, the solver becomes the space-marching scheme. 
Under this modified scheme, a solution is obtained plane 
by plane from upstream to downstream by carrying out a 
sufficient number of implicit iterations in each plane 
until the convergence criterion is met. A solution for the 
entire computational domain is established in a single 
forward sweep. 

The t-e Turbulence Model 

The governing equations of the Jones-Launder for- 
mulation of the k-£ turbulence model are written as a pair 
of coupled transport equations in conservative form. In 
principle, this pair can be implemented together with the 
Navier-Stokes equations as either a set of seven coupled 
equations or a separate pair uncoupled from the Navier- 
Stokes equations. The fully coupled approach would 
cause serious problems such as a significant increase in 
computational effort and working space in the computer 
memory and numerical stiffness of the coupled set of 
seven equations. In the PAB3D code, solutions of the k 
and e equations are decoupled from the Navier-Stokes 
equations and from each other. Time step differences 
remain in this uncoupled system of flow and turbulence 
equations. However, the problem is circumvented by 
solving these k and e equations with a CFL (Courant, 
Fredricks, and Levy) number that is smaller by at least 
a factor of 2. The potential difference in timewise devel- 
opment of the flow variables and turbulence variables 
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has not led to any noticeable effect in either the overall 
convergence rate or the quality of the solutions. 

Multiblock Structure and Boundary Conditions 

The PAB3D code is designed to handle complex 
configurations by using several types of multizone, 
multiblock grid topologies. A restricted option, which is 
suitable for jet plume calculation with the space- 
marching schemes, calls for streamwise division of the 
computational domain into zones in the ^-direction. The 
grid space in each zone can further be divided into blocks 
in the r|- and ^-directions. Otherwise, the code supports a 
general multiblock scheme where the computational 
domain can be divided into any collection of blocks. 
Number of blocks, block size, and parametric orientation 
are not restricted. The concept of zones is not relevant in 
this general scheme. General patched block interface 
communication is allowed. The only restriction for this 
general multiblock connectivity scheme is that the con- 
nected block interfaces are contiguous. A grid partition 
feature is available in the PAB3D code for the conve- 
nience of turbulence modeling. If different viscous stress 
models are employed within a block, the ^-direction of 
the block can be partitioned by choosing a starting index 
for each viscous stress model. 

The boundary conditions often used for jet computa- 
tions include inflow, outflow, free stream, solid walls, 
and geometrical symmetry. Three types of inflow bound- 
ary conditions are provided: Riemannian characteristics, 
fixed total temperature and total pressure, and a com- 
pletely fixed set of five flow parameters. Two outflow 
boundary conditions are needed: constant pressure for 
subsonic flows and extrapolation for supersonic flows. 
The Riemannian characteristics boundary condition is 
used at free stream boundaries. On a solid boundary, 
either a no-slip or an inviscid-slip boundary condition 
can be specified. Finally, the symmetry boundary condi- 
tions include mirror imaging across a plane in any orien- 
tation and polar symmetry around an axis in the 
streamwise direction. 

In addition, a universal high-order symmetry bound- 
ary condition for Navier-Stokes code applications is 
developed in the course of this jet plume study. This uni- 
versal symmetry boundary condition provides a simple 
method for the user to specify a symmetry boundary con- 
dition at a grid plane not aligned to a surface with a con- 
stant physical coordinate value. Reference 14 details this 
boundary condition. 

Adaptive Grid Algorithm in the PAB3D Code 

For nonaxisymmetric jet calculations in this report, 
one quarter of the jet cross section is represented in the 
grid domain. Hence, flow symmetry across both the hori- 


zontal and the vertical axis is assumed. In each plane, the 
grid is divided into two parts: a high density grid in the 
near field of the jet flow and a low density grid for the far 
field. Only the high density grid near the jet flow is 
adapted to the flow solution. A Cartesian topology is 
chosen for the initial unadapted high density grid. Grid 
adaption proceeds from plane to plane in the streamwise 
direction. Two monitoring functions are used for adapt- 
ing the grid to the velocity and pressure gradients of the 
flow solution. The monitoring functions are normalized 
so that one constant is used for each function to control 
the intensity of adaption. (See ref. 21 .) A single grid was 
used in reference 21 to cover both the near field and the 
far field. A third monitoring function was employed to 
redistribute a uniform Cartesian grid to form a dense grid 
zone in the near field and a sparse grid distribution in the 
far field. 

The adaptive grid procedure is coupled to the space- 
marching solver in PAB3D. Grid indices i, j, and k are 
assigned to the T|, and £ coordinates of the generalized 
coordinates. In the space-marching algorithm, a numeri- 
cally converged solution is computed at each j % /c-plane 
through multiple iterations. This solution is then coupled 
to the next plane downstream, and the computational 
process is repeated. 

For jet flow computations considered in this report, 
initial conditions for the first plane representing the flow 
condition at the nozzle exit are prescribed according to 
known nozzle operating conditions. The initial grid at the 
first plane is generated externally according to the initial 
flow conditions by using the same grid adaption proce- 
dure as the one implemented in the PAB3D code. Once 
the solution process is started, grid adaptions for subse- 
quent grid planes are computed within the code. A grid is 
created for the (/ + 1 )-plane by adapting the grid to the 
numerically converged solution in the /-plane. The adap- 
tive grid procedure is implemented as efficiently 
as possible to match the high efficiency of the space- 
marching solver. The computational efficiency of this 
multiple function grid adaption procedure is analyzed 
during this study. The time taken for grid adaption is 
approximately 4 percent of the total time required for the 
flow solver. Only 1 cycle of adaption is used for each 
plane. The flow solution at each plane normally takes 
approximately 20 to 30 iterations before nominal conver- 
gence criterion is met. 

Results and Discussion 

Results of jet flow computations are presented 
in several groups: jets operating at on-design conditions, 
supersonic jets containing weak shocks, supersonic jets 
containing strong shocks, and nonaxisymmetric super- 
sonic jets. For these groups of computation, the initial 
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flow condition for the jet is specified at the inflow 
boundary of the computational grid. The internal flow 
upstream of the jet nozzle is not modeled. Figure 1 shows 
a sketch of the jet flow configuration and a typical com- 
putational grid for on-design axisymmetrical jets. 

On-design operation of a jet is defined as the condi- 
tion for which the jet exit static pressure is identical to 
the ambient static pressure. For a subsonic jet, the exit 
static pressure is naturally adjusted to the ambient static 
pressure. A supersonic jet flow is established by using a 
convergent-divergent nozzle designed for a fixed Mach 
number. The on-design nozzle pressure ratio (NPR), 
defined as pjp 0 , is a fixed value for each given Mach 
number. For jets operating at on-design NPR, pressure 
gradients are very small in the entire flow domain. 
Shocks in the flow domain are typically weak or absent 
in on-design supersonic jets. The principal driving mech- 
anism for on-design jet plume development is turbulent 
mixing in the jet shear layer. This report computes on- 
design jet flows by using the space-marching scheme in 
PAB3D. 

At off-design operating conditions, the initial jet 
flow condition is either overexpanded or underexpanded. 
Shock waves will appear in the jet flow. For a given noz- 
zle geometry, the exit jet Mach number is fixed regard- 
less of NPR (assuming that the NPR is high enough to 
fully establish supersonic flow at the nozzle exit). At 
NPR values sufficiently close to the design point, only 
weak shocks are present in the jet flow and jet flow 
development can be computed using the space-marching 
solver in PAB3D. 

Once strong shocks in the form of Mach disks appear 
in the jet plume, the flow downstream of the shock 
becomes subsonic. Furthermore, high static pressure 
immediately downstream of the shock leads to rapid 
acceleration and expansion of the subsonic flow. Hence, 
a strong pressure gradient exists. The time-marching 
solvers in the PAB3D code must then be used because 
conditions permitting the use of the space-marching 
scheme are violated in this region. However, the space- 
marching method alone cannot detect the occurrence of a 
Mach disk. The decision to use either the space-marching 
or time-marching options in the PAB3D code for a par- 
ticular case must be guided by tabulated experimental 
data or by theoretical estimates. Reference 5 gives an 
excellent reference for Mach disk formation in axisym- 
metric jet plumes. 

Jet exhaust nozzles of practical interest in propulsion 
systems may have a nonaxisymmetric exit cross section. 
The dynamic characteristics of nonaxisymmetric jets are 
significantly more complex than those for axisymmetric 
jets because of the added degrees of freedom in jet flow 
geometry. The adaptive grid option in the PAB3D code 


is used to provide appropriate grid density distribution 
for the shear layer and shock regions in the jet flow. The 
following subsections give detailed discussions of results 
in each group of jet flow computations. 

On-Design Circular Jet Plumes 

Flow solutions for on-design circular jets within jet 
exit Mach numbers ranging from 0.6 to 3.0 are computed 
using the space-marching method in the PAB3D code. In 
this series of jet flow simulations, the initial flow condi- 
tion for the jet plume is specified at the inflow boundary 
of the computational domain. A small velocity compo- 
nent in the ambient air parallel to the jet flow direction is 
required in the PAB3D code for maintaining numerical 
stability of the space-marching scheme. A freestream 
Mach number of 0.001 is sufficient in fulfilling this 
numerical requirement. 

The on-design jet grid is constructed as a single layer 
wedge which covers a sector of 2.5° in the circumferen- 
tial direction. Figure 1 shows general layout of this grid. 
There are 400 uniformly sized grid cells in the stream- 
wise direction covering a distance of x/R = 40 and 48 
grid cells in the transverse direction covering a radial dis- 
tance of y/R = 8. At the inflow station of the jet, the jet 
plume is defined by 18 grid cells, and the remaining 30 
grid cells cover the distance from y/R = 1.0 to 8.0. The 
initial shear layer region near the nozzle exit plane is 
covered by 24 grid cells centered above and below the 
nozzle lip. High grid density is provided in the shear 
layer to capture the turbulent mixing process. As the jet 
flow spreads downstream, approximately 30 grid cells 
are located within the jet flow. For computational conve- 
nience, the grid domain is divided into four blocks in the 
streamwise direction. The grid domain can easily be 
extended in the streamwise direction by adding more 
blocks. 

General features of jet flow computation using the 
PAB3D code with the two-equation k-e turbulence clo- 
sure model are illustrated by the solutions of a typical 
subsonic jet at M — 0.6. Compressibility correction for 
the k-e turbulence model is not needed in this com- 
putation. Figure 2 shows the computed centerline veloc- 
ity profile for the M - 0.6 jet. The centerline velocity 
maintains its exit value for a distance up to approxi- 
mately x/R =12 and decreases farther downstream as a 
result of turbulence mixing. The classical relation of 
velocity decay is given by 


where L c is the intercept of the x 1 decay curve and hori- 
zontal line u ( (x)/U € = 1 .0 (referred to as the potential 
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core length in jet flow literature). Figure 2 shows both 
the computed centerline velocity profile and the classical 
velocity decay as indicated by equation (6). Good agree- 
ment is shown between the PAB3D solution and the 
result obtained with equation (6). 

Figure 2 also shows the computed centerline velocity 
profiles for a Mach 2.0 jet operating on-design using the 
Jones-Launder k-E turbulence model with three different 
methods of compressibility corrections. The compress- 
ibility correction factor in the k-E turbulence model has a 
strong influence on jet flow development. Turbulence 
mixing is strongest in the jet flow when no compressibil- 
ity correction is applied. For this case, the potential core 
length is LJR = 17.2. The action of turbulence mixing in 
the jet is weaker when compressibility corrections are 
applied. The value of LJR is 22.6 and 25.2 for the Sarkar 
and the Wilcox methods, respectively. The velocity 
decay downstream of the end of the potential core is also 
computed according to equation (6) and the value of LJR 
for each case. The results are shown in figure 2. 

Good agreement is observed between the PAB3D 
solutions using the Sarkar and the Wilcox corrections 
and their corresponding results using equation (6) for 
uJU e greater than 0.7. For uJU e less than 0.7, the 
PAB3D solution begins to deviate from the classical \/x 
decay rate. For the solution without compressibility cor- 
rection, the decay rate starts to deviate from the \/x decay 
at approximately uJU € - 0.8. Without compressibility 
correction in the k-E turbulence model, the predicted tur- 
bulence level is too high for the Mach 2.0 supersonic jet 
solution. This steep velocity decay is an indication of 
excess mixing in the jet shear layer. 

Figure 3 shows the downstream evolution of the 
M = 0.6 jet velocity cross section. At xIR = 0, the initial 
velocity profile across the entire width of the jet nozzle 
exit has a prescribed constant value of U e . The cross sec- 
tion at x/R = 5 (fig. 3) shows the initial development of a 
thin shear layer, and the width of the potential core is 
narrower than its width at the jet exit. The cross section 
at x/R = 15 is located just downstream of the end of the 
potential core. The velocity profile at x/R = 15 has not 
yet attained a Gaussian distribution. However, the Gauss- 
ian velocity distribution has been established at x/R = 25. 
Figure 4 shows the turbulence intensity distributions at 
the corresponding x/R stations. The peaks of the turbu- 
lence intensity distributions at x/R = 5 and 15 are located 
in the middle of the shear layer where the velocity gradi- 
ent is the highest. Although the centerline turbulence 
level at x/R = 25 is significantly higher, the peak turbu- 
lence intensity remains off center, and the turbulence 
intensity distribution across the jet is not Gaussian. 
According to measured data by Wygnanski and Fiedler 


(ref. 22), self similarity of the turbulence intensity is usu- 
ally established at x/R values between 50 and 70. 

Figure 5 shows computed centerline velocity profiles 
for a M = 2.22 jet and the experimental data measured by 
Eggers. (See ref. 23.) Like the centerline velocity profiles 
for a Mach 2.0 jet shown in figure 2, the solutions 
obtained by using different compressibility corrections 
are different. With no compressibility correction, the 
potential core length is underpredicted. The location of 
the end of the potential core appears to agree with the 
centerline velocity profile predicted using the Wilcox 
model. However, centerline velocity computed by using 
equation (6) and the potential core length of the Sarkar 
solution L c /R = 27.15 agrees very well with the data 
obtained farther downstream. (See ref. 23.) The agree- 
ment between computational and measured data is much 
better when compressibility corrections are applied, 
although a small difference exists between the Sarkar 
model and the Wilcox solutions. Figure 6 shows the cor- 
responding results of velocity distributions in the jet 
cross section at x/R = 25. The importance of 
compressibility correction for supersonic jets is further 
illustrated here, as the compressibility-corrected compu- 
tations come very close to the measured data, whereas 
the uncorrected computation underpredicts the centerline 
velocity by nearly 40 percent. 

A group of on -design jet plumes with exit Mach 
numbers ranging from 0.6 to 3.0 is computed to illustrate 
the trend of turbulent mixing as a function of Mach num- 
ber. Figure 7 shows typical turbulence intensity distribu- 
tions u/U e in the longitudinal plane of symmetry of the 
jet at three different Mach numbers: 0.8, 1.2, and 1.6. 
The contours in figure 7 show that turbulence is absent in 
the potential core region. Intense levels of turbulence 
start to develop at the lip of the jet nozzle exit. The posi- 
tion of the maximum turbulence intensity in the initial 
zone of the shear layer occurs near the lip line of the jet. 
As the shear layer evolves farther downstream, the posi- 
tion of maximum turbulence intensity migrates towards 
the jet centerline. This general pattern remains the same 
for all on-design circular jets computed within the Mach 
number range from 0.6 to 3.0. The length of the potential 
core and the value of maximum turbulence intensity vary 
as a function of Mach number. Figure 8 summarizes the 
computed turbulence intensity distributions along the jet 
centerline as a function of Mach number. The Wilcox 
compressibility correction is used for these computations 
because the definition of the Wilcox correction provides 
a consistent blending of compressibility correction for 
subsonic and supersonic flow regions. 

The potential core length is usually defined as the 
distance from the jet exit to the beginning of centerline 
velocity decay. An important equation for potential core 
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length as a function of jet Mach number is given by Lau, 
Morris, and Fisher in reference 24 as: 

L c 

-g = 8.4 + 2.2M 2 (7) 

Core length is obtained by an empirical curve fit to a 
large collection of measured values for potential core 
length in subsonic and supersonic jets. 

Figure 9 shows the potential core lengths computed 
with the standard Jones-Launder two-equation k-E turbu- 
lence model with the Sarkar and the Wilcox compress- 
ibility corrections. The core length is defined as the point 
where the value of the centerline jet velocity has dropped 
to 0.99 times the jet exit velocity. The potential core 
length derived with the Wilcox compressibility correc- 
tion is higher than the value computed in equation (7) for 
the entire Mach number range. Values obtained by using 
the Sarkar compressibility correction are higher than the 
Wilcox results. However, the trends of core length varia- 
tion as a function of Mach number are similar in all three 
sets of results. An alternate value of the potential core 
lengths can be obtained from the computed jet flow solu- 
tion when the end of potential core in the jet flow is 
defined as the point where turbulence intensity level 
exceeds a threshold of = 0.01 along the jet center- 
line. The core lengths obtained by the turbulence inten- 
sity definition (also shown in fig. 9) agree very well with 
the values obtained with the velocity decay criterion. 

The difference between the computed and the empir- 
ical curve fit formula based on measured values origi- 
nates from several sources. In the work by Lau, Morris, 
and Fisher, the experimental database contains measured 
potential core length values for jets operating at different 
temperatures. Equation (7) is a curve fit for isothermal 
jets where the jet static temperature is the same as the 
ambient temperature, whereas the jet total temperature is 
higher than the ambient total temperature. Many data 
points for cold jets, where the jet total temperature is the 
same as the ambient air temperature and therefore the jet 
static temperature is colder than the ambient temperature, 
are above the curve fit of equation (7). In this report, the 
jet flows are computed as cold jets. 

A second source of discrepancy may come from 
modeling boundary conditions in the computations. For 
jet flows in the laboratory, the boundary layer within the 
jet nozzle has a finite thickness at the nozzle exit. The 
initial turbulence level and the thickness of the nozzle 
boundary layer give the jet mixing layer an earlier start in 
its development. Therefore, the computed core length 
will be shorter if the initial boundary layer at the jet noz- 
zle exit is included in the computations. In addition to 
these circumstances, grid density and accuracy of turbu- 
lence modeling are important factors to be considered for 


further refinements of the computational method for jet 
flow predictions. 

It is significant that the mean flow and turbulence 
levels of on-design circular jet plumes are predicted over 
such a wide range of Mach numbers by using the stan- 
dard Jones-Launder &-£ turbulence model and the Wilcox 
compressibility correction without changing the pub- 
lished constants for the turbulence model. In a broader 
context, the modeling of jet plumes is often required in 
propulsion and airframe integration. A consistent compu- 
tational analysis for such jet plume modeling using the 
Navier-Stokes method should not permit ad hoc changes 
to the turbulence model. The results of this parametric 
investigation indicate that ad hoc modifications to the 
standard Jones-Launder turbulence model are not needed 
for jet flow analysis. 

Off-Design Jets Containing Weak Shocks 

This section shows the flow properties of a Mach 2 
jet containing weak shocks using solutions obtained 
within a limited range of nozzle pressure ratios. The 
space-marching solver procedure in the PAB3D code is 
used to compute these jet flows. At Mach 2, the jet flow 
is free from Mach disk formation for values of NPR 
between 4.6 and 13.8, which correspond to a ratio of jet 
exit to ambient pressure p e /p Q between 0.6 and 1 .8. Fig- 
ure 10 shows a density contour for a typical under- 
expanded jet. At the jet exit, a curved shock near the lip 
line of the jet nozzle is formed to resolve the pressure 
difference between the ambient flow and the flow inside 
the jet. An internal weak shock system which reflects 
repeatedly between the shear layer and the jet centerline 
also exists in the jet. 

Figure 1 1 shows the computed pressure distribution 
along the jet axis for p e /p 0 = 0.8, 1.2, 1.4, 1.6, and 1.8. 
Only one overexpanded case is included in this collection 
(P e /Po = ®‘ 8)« The pressure distribution of the over- 
expanded jet is characterized by a sharp shock front and 
pressure peak close to the jet exit. This feature is not 
found in underexpanded jets. Beginning with the second 
peak, the features of overexpanded jet pressure oscilla- 
tions on the jet centerline follow the same trend as the 
patterns shown for underexpanded jets. For the under- 
expanded jets, the cycle of pressure oscillation begins 
with a smooth expansion. The flow is then recompressed 
towards the first pressure peak in two stages. The pres- 
sure rises sharply about halfway then recompresses grad- 
ually the rest of the way. Although not shown in 
the figure, expansion and recompression processes are 
smooth for subsequent periods of oscillation. Refer- 
ence 25 gives further discussion of the centerline pres- 
sure distributions of off-design jet flows at Mach 2. 


10 



Extensive flow visualization measurements of super- 
sonic jets at off-design conditions were obtained by Love 
et al. (See ref. 5.) Shock formation in the jet flow is char- 
acterized by two lengths: l h the distance from the jet exit 
to the first shock intersection with the jet centerline; and 
W , the distance from the jet exit to the first shock inter- 
section with the sonic line in the shear layer. Figure 10 
shows definitions of these lengths. Figure 12 shows the 
computed values and measured values of l\ and W 
(ref. 5) for several overexpanded and underexpanded val- 
ues of NPR. Excellent agreement is demonstrated by the 
results in figure 12. 

Figure 13 details computed centerline pressure dis- 
tributions and experimental data for a Mach 2.0 jet 
at p e !p 0 = 1 -445 (NPR =11.3). Three solutions are 
obtained by using the basic k-E turbulence model with no 
compressibility correction, the Sarkar correction, and the 
Wilcox correction. The solution without compressibility 
correction shows that the amplitude of pressure oscilla- 
tions diminishes rapidly downstream and predicted 
wavelength is much shorter than the experimental data in 
the downstream region of the jet. The solutions obtained 
with a compressibility-corrected k-E turbulence model 
show excellent agreement with measured data. Differ- 
ences between the Sarkar and Wilcox corrections are 
small. The amplitudes of the computed solutions closely 
follow the test data, but their phase relations with re- 
spect to the measured data are somewhat different. 
At xIR = 40, the Wilcox solution leads the measured data 
by approximately one sixth of one period, whereas the 
Sarkar solution lags behind the measurements by approx- 
imately half that amount. All three solutions are very 
similar near the jet exit. However, the amplitude of the 
first pressure peak near the jet exit is underpredicted by 
approximately 15 percent. 

Figure 14 shows the computed values for the axial 
turbulence velocity component and the measured data 
obtained by Seiner, Dash, and Wolf (ref. 26), The inter- 
action between the repeated shock-cell structure and the 
turbulence produces a significant periodic modulation of 
the axial turbulence velocity component. The magnitude 
of the fluctuation is in phase with the pressure fluctuation 
in the jet. (See fig. 13.) Good agreement in both phase 
and amplitude is seen between the computed solutions 
and the measured data. The compressibility-corrected 
solutions provide better agreement with the measure- 
ments than the uncorrected solution. It is encouraging to 
find from this comparison that the standard Jones-Laun- 
der k-E model is capable of accurate predictions of the 
turbulence velocity in a shock-containing supersonic jet. 
For practical applications such as jet noise prediction, an 
estimate of turbulence intensity in the jet flow is needed. 
A computational capability for predicting turbulence 
intensity distributions in the jet flow is highly desirable 


because measurement of turbulence in high speed flow is 
exceedingly difficult. 

Better predictions of the turbulent velocity fluctua- 
tions in a supersonic jet can be obtained with further 
improvement of the turbulence model. In standard k-E 
turbulence models, the local turbulence kinetic energy is 
attributed equally to all three turbulence velocity compo- 
nents. However, it is known that the magnitude of the 
axial component is higher than those of the transverse 
components in the jet shear layer. Therefore, a better 
redistribution relationship of the turbulent kinetic energy 
and the Reynolds stress tensor components would raise 
the value of the computed uJU £ . Furthermore, the mod- 
ulation of the axial turbulence velocity component by the 
internal shock waves would be stronger, since the ampli- 
tude of shock turbulence interaction is roughly propor- 
tional to the shock strength and the magnitude of the 
axial component of the velocity fluctuations. 

Use of the space-marching algorithm to obtain a jet 
flow solution requires less than 100 seconds of CPU time 
on the Cray Y-MP computer at the Langley Research 
Center. Use of the time-marching solver to obtain a con- 
verged solution for the same cases typically requires 
2000 seconds of CPU time. The ratio of computer time 
required when using the time-marching solver increases 
by a factor of 20. Figure 15 presents jet centerline pres- 
sure distributions obtained by using the space-marching 
and time-marching solvers. The flow solutions obtained 
by these two different procedures are practically the 
same. Detailed discussion of this comparison can be 
found in reference 27. 

Off-Design Supersonic Jets Containing a Mach 

Disk 

Mach disks may appear in a supersonic jet if the 
operating NPR is significantly different from NPR 
design value. The conditions for Mach disk formation for 
a given nozzle depend on nozzle design Mach number 
and details of the nozzle geometry, such as the nozzle 
wall exit angle. Mach disk formation can occur in both 
overexpanded and underexpanded conditions. For a 
Mach 2.0 nozzle with on-design NPR of 7.82, Mach disk 
appears if the operating NPR is less than 4.6 or greater 
than 13.8. For a nozzle with M = 1.5 with an on-design 
NPR value of 3.67, a Mach disk will form in the jet for 
NPR less than 2.7 or greater than 6. 1 . 

For jet flow computations where Mach disk forma- 
tion is expected, the time-marching solver in the PAB3D 
code is used. A different computational grid is also 
required. When the case of a Mach disk containing jet 
flow is originally computed with the on-design jet grid, 
the Mach disk is never formed in the converged solution. 
The shock front initiated at the nozzle exit propagates as 
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a weak shock all the way to the jet centerline and then 
continues as a regular reflected weak shock. In the 
on-design grid, the cell streamwise versus radial aspect 
ratio is 4. Although the PAB3D code solver is designed 
as an upwind algorithm, certain numerical errors in the 
transonic regime prevent the proper formation of a Mach 
disk in the flow solution. Because general patched grid 
capability is available in the PAB3D code, a new grid is 
created so that the grid for x/R from 0 to 2 has a grid dis- 
tribution similar to the on-design jet grid but with double 
the density in each direction. For x/R from 2 to 10, a 
uniformly sized grid distribution is retained in the 
streamwise direction. In the radial direction, a uniform 
grid distribution is provided for y/R from 0 to 2 so that 
cell aspect ratio in the entire region is 1 .0. An exponen- 
tially expanding grid is used for y/R from 2 to 8 to cover 
the free-stream domain outside the jet flow. Figure 16 
shows a sketch of this revised grid. The overall grid 
domain is divided into four blocks for computational 
convenience. The correct Mach disk containing solution 
is obtained by using this revised grid. 

A solution for a Mach 1 .5 jet operating at 
NPR =11.6 (p e /p 0 = 3.15) is obtained. Figure 17 shows 
the density and Mach number distributions in the jet. A 
well-formed Mach disk is located at x/R = 4.4. The 
radius of the Mach disk is approximately 0.68/? in the 
computed solution. The location of this Mach disk agrees 
with the measurements given in reference 5. However, 
the computed radius of the Mach disk is smaller than the 
corresponding measured value. The reflected weak shock 
and a slip line initiated at the outer edge of the Mach disk 
is clearly shown by the computed density contours. The 
contour value indicates that the Mach number upstream 
of the Mach disk has accelerated to values greater 
than 4.0, whereas the Mach number downstream of the 
first Mach disk is reduced to values below 0.2. Down- 
stream of the first Mach disk, the flow near the centerline 
again accelerates to supersonic speeds near x/R = 8.0. A 
second Mach disk is subsequently formed at x/R = 8.6. 
Though much weaker, the second Mach disk can be seen 
in a schlieren photograph for a jet operating at nearly the 
same jet initial conditions. (See ref. 5.) 

Since the time-marching computations for the Mach 
disk case are executed by using grid sequencing, con- 
verged solutions at three grid levels are obtained during 
the process. Figure 1 8 shows the Mach number contours 
using the one-fourth and one-half linear grid density in 
the y- and k-directions. Even at the quarter density grid 
level, the first Mach disk is captured in the solution. Both 
the location and the radius of this Mach disk are estab- 
lished in this coarse grid. In the half density grid, the sec- 
ond Mach disk emerges in the solution. Only minor 
changes in the flow physics are detectable between the 


half density grid solution and the full density grid solu- 
tion, as shown earlier in figure 17. 

Adaptive Grid Computations of 
Nonaxisymmetric Jets 

For jet flow computations using an adaptive grid, a 
quarter plane symmetry for the jet is assumed. The grid is 
divided into two domains: a high grid density inner 
domain near the jet flow, and an outer domain with 
reduced grid density to cover the free-stream domain 
away from the jet flow. A Cartesian grid topology is used 
in the inner domain to accommodate a wide range of jet 
exit geometries. The outer domain has a polar topology 
with significantly less grid density than the inner domain. 
Figure 19 shows a sketch of the grid cross section. Refer- 
ence 28 shows that the computational simulation of a cir- 
cular jet remains perfectly axisymmetric even though the 
grid is Cartesian. Furthermore, the adaptive grid proce- 
dure provides adequate grid densities to support accurate 
computations in the jet shear layer and in regions near a 
shock front. 

This section discusses computed solutions for an 
elliptic, a rectangular, and a square jet using the adaptive 
grid method. The elliptic jet is known for its unusual 
mixing characteristics. The rectangular jet family, which 
includes the square jet as a special case, is widely used 
for propulsion integration in advanced aircraft systems. 
Both elliptic and rectangular jets are capable of switching 
their major and minor axis directions in different cross 
sections along the jet. An initial Mach number of 2.0 and 
an operating NPR of 1 1.31 are chosen for all three con- 
figurations. The operating NPR corresponds to an exit 
static pressure ratio of p e /p 0 = 1.445. In addition, the 
Mach 2 elliptic jet is also computed at its on-design flow 
condition at NPR = 7.82. 

Figure 20 shows the computed Mach number con- 
tours for the underexpanded elliptic jet in the major and 
the minor planes of symmetry. The elliptic cross section 
at the nozzle exit plane has an aspect ratio of 2.0. The 
shock fronts are clearly defined by the Mach number 
contours in the major plane of symmetry. The shock 
reflection pattern in the core region of the jet is quite reg- 
ular. In the plane containing the minor axis, the shock is 
initially reflected in the shear layer with a scale propor- 
tional to the minor axis length. However, the short wave 
pattern quickly disappears at approximately x/R = 6. Far- 
ther downstream in the jet, only the long wave pattern 
dominated by the length scale of the major axis remains. 
Another unique feature in figure 20 is the expansion rate 
of the outer jet boundary. In the plane containing the 
major axis, the jet boundary expands slowly in the radial 
direction. In contrast, the jet boundary in the plane 
containing the minor axis expands rapidly in the radial 
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direction. By approximately x/R = 15, the width of the 
elliptic jet in the original minor axis direction has 
exceeded the jet width in the original major axis direc- 
tion. Hence, this computation indicates an axis switching 
phenomenon for a supersonic elliptic jet. 

When using similarity analysis, the major and the 
minor axes of the initial cross section are considered two 
independent reference scale lengths for the jet flow. A 
simple consequence of assuming two independent refer- 
ence scales would be that the internal shock reflection 
pattern would repeat in two directions along two differ- 
ent scales. On the other hand, the difference between 
these two scales must be resolved within the jet flow to a 
common scale since shock fronts cannot cross each other 
without some type of interaction. The computed result 
demonstrates the complexity of such aerodynamic inter- 
actions. In the plane of symmetry containing the minor 
axis, the internal shock wave length is initially governed 
by the minor axis scale. However, the shear layer posi- 
tion expands rapidly outward in the minor plane of sym- 
metry; thus, the reflection length scale of the downstream 
shock wave pattern is changed. In the major plane of 
symmetry, the width of the jet in this plane remains 
approximately constant; thus, the reference scale of the 
jet in the minor axis direction is allowed to catch up. The 
nonlinear interaction within the jet flow eventually leads 
to a unified scale length for the shock cell system. 

Figure 21 shows the Mach number contours in an 
on-design elliptic jet at Mach 2.0. In the absence of a 
shock structure in the jet, the Mach number distribution 
in both the major and minor planes of symmetry is 
smooth and indistinguishable from the previously com- 
puted Mach number distributions in circular jets. Similar 
to the underexpanded elliptic jet case, the shear layer 
growth in the minor plane of symmetry is faster than the 
growth in the major plane of symmetry. At x/R = 40, the 
widths of the jet as seen in both planes of symmetry are 
almost equal. However, axis switching does not occur in 
the on-design case. 

In order to examine the possibility of axis switching, 
cross section Mach number contours are shown for these 
four jets in figures 22-25. The exit cross section aspect 
ratio for the elliptic and the rectangular jet exit shapes 
is 2.0. The Mach number contours in each of the 
cross sections show only a narrow band from M = 0.8 
to M = 1.2 with a contour interval of 0.1 to highlight the 
shape of the cross section. Figure 22 shows the evolution 
of the elliptic underexpanded jet cross sections. At the 
exit, the major axis of the elliptic cross section is oriented 
in the horizontal direction. The jet boundary grows 
rapidly in the vertical direction. At x/R = 30, the major 
axis of the elliptic cross section has clearly switched 


to the vertical direction. The aspect ratio of the ellipse at 
x/R = 30 is approximately 1.50. 

Figure 23 shows the equivalent sequence for a rect- 
angular jet. Jet boundary growth in the vertical direction 
is even faster than that of the elliptic jet. At x/R = 30, the 
aspect ratio of the shape is approximately 1 .70. The cor- 
ners of the original rectangular shape have been rounded 
off in the process of turbulent mixing. 

Figure 24 shows the evolution sequence for a square 
jet. In this case, the original square shape for the Mach 
number contours evolves rapidly in the jet flow. At 
x//?= 12, the corners of the square are actually trans- 
posed by 45°, with the comer sharpness well preserved. 
As the shear layer grows thicker farther downstream, the 
shape of the jet cross section quickly losses its distinction 
as a square and eventually becomes circular. 

Figure 25 shows the shape evolution sequence for an 
elliptic jet with on-design exit pressure ratio. Although 
the jet grows mainly in the vertical direction, axis switch- 
ing does not occur in this on-design jet flow. The jet sim- 
ply becomes a near-circular jet at x/R = 30. At the last 
computed streamwise position at x/R = 40, the cross sec- 
tions for all four jets simply retain their geometrical char- 
acters similar to those established at x/R = 30. 

Figure 26 shows the Mach number contours in the 
plane of symmetry of the square jet and its centerline 
Mach number distribution. Qualitatively, these distribu- 
tions are very similar to the corresponding distributions 
shown in figure 27 for a circular jet operating at the same 
exit Mach number and NPR. The visually striking 
dynamic behavior of the square jet axis (fig. 24) shows 
that switching apparently has little influence on flow 
development near the jet centerline. It should be pointed 
out that the square jet and the circular jet, with a common 
reference dimension of 1.0, have different jet exit areas. 
In order to compare the streamwise jet flow development 
on a normalized basis, the length scale for the square jet 
in figure 24 should be reduced by a factor of 

Jn/A = 0.8862. With this scale adjustment, the phase 
and amplitude of the centerline Mach number oscilla- 
tions of the square jet and the circular jet agree almost 
exactly starting from the second peak. 

In order to provide some validation for the adaptive 
grid computation procedure for nonaxisymmetric jets, a 
Mach 2.0 on-design circular jet solution computed by 
using the axisymmetric grid is compared with the same 
jet computed by using a three-dimensional adaptive grid. 
Two levels of adaptive grid densities are also used to ver- 
ify grid convergence: 40 x 40 cells and 56 x 56 cells for 
the inner high density grid cross sections. Figure 28 
shows the axisymmetric grid solution for M = 2.0. The 
adaptive grid results are shown in figures 29 and 30. 
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Nearly identical solutions for the turbulence intensity 
distribution in the meridian plane and along the center- 
line of the jet are obtained for the two adaptive grid 
densities. For example, the maximum turbulence level in 
the jet plume is 0.134 for both the 40 x 40 cell and the 
56 x 56 cell solutions. The solution using an axisymmet- 
ric grid differs slightly from the adaptive grid results in 
two aspects. First, the maximum turbulence level in the 
middle of the shear layer is slightly higher, with a value 
of 0.142. Second, the turbulence intensity profile along 
the centerline is shifted upstream by approximately 
xIR = 2.0. 

The reason for the different maximum intensity in 
the shear layer is not clear. However, the spatial shift of 
the turbulence intensity profile along the centerline has a 
geometrical explanation. It is difficult for the grid adap- 
tion algorithm to handle very large velocity gradients 
such as those occurring near the jet exit. Consequently, it 
is not possible to specify an initial shear layer thickness 
of less than 0.05 jet radius at the jet exit plane. Since the 
initial shear layer is thicker, the inner boundary of the 
turbulent shear layer intersects with the centerline at a 
smaller value of x/R. In fact, with a downstream shift of 
the adaptive grid centerline turbulence intensity profile, 
the turbulence intensity profile can be matched perfectly 
with results using the axisymmetric grid. 

Concluding Remarks 

The main purpose of this report is to establish a uni- 
fied method for jet flow prediction using the Navier- 
Stokes method with a two-equation k-z turbulence clo- 
sure model. Although the jet flow may contain a variety 
of complex flow physics features, the Navier-Stokes 
method simply requires that the initial condition and 
boundary conditions of the jet operating conditions be 
specified for the problem. Detailed flow physics devel- 
opments in the jet are predicted by the Navier-Stokes 
method. The validity of this approach is demonstrated by 
the high quality jet flow solutions obtained with the 
PAB3D code. 

This study examines several categories of jet flow 
conditions. For on-design subsonic and supersonic axi- 
symmetric jets, the flow field is dominated entirely by 
turbulent mixing. Numerical solutions within a Mach 
number range of 0.6 to 3.0 are accurate when compared 
with available experiment data for parameters such as 
mean velocity and turbulence intensity distributions in 
the jet, centerline velocity decay, and the potential core 
length variation as a function of Mach number. 


For off-design supersonic jet flows containing weak 
shocks, flow predictions are compared with experimental 
data. Good agreement is obtained between the computed 
results and experimental data for key parameters, includ- 
ing first shock-cell lengths and centerline pressure distri- 
bution. The predicted distributions of the streamwise 
component of turbulence velocity fluctuation in an 
underexpanded Mach 2.0 jet show good agreement with 
measured data. 

Turbulence intensity in the jet flow is an important 
quantity for jet noise prediction. Since direct measure- 
ment of turbulence in a supersonic jet is very difficult to 
make, a predictive capability provided by the PAB3D 
code is very useful for practical applications. Good 
agreement between predictions and experimental mea- 
surements has also been obtained for a Mach 1 .5 jet oper- 
ating at 3.15 times its design nozzle pressure ratio where 
Mach disks are present in the jet flow. 

Many of the modem propulsion jet nozzles employ 
nonaxisymmetrical exit geometries. The adaptive grid 
method examined in this study has produced good results 
for elliptic, rectangular, and square jets. However, the 
computed results are not verified for lack of experimental 
data. The accuracy of the adaptive grid procedure is illus- 
trated by a comparison between an adaptive grid solution 
of an axisymmetric jet and a solution for the same jet 
using a single cell wedge grid. Although the adaptive 
grid has a Cartesian-topology and the single-cell wedge 
grid has cylindrical symmetry boundary conditions, the 
solutions are essentially identical. 

For most jet flows where strong shocks are absent in 
the computational domain, the space-marching solver in 
the PAB3D code can be used. When the space- marching 
option is used for jet flow computation as conditions per- 
mit, the computer time is one twentieth of the time 
required for obtaining a time-marching solution with the 
same flow conditions. The accuracy of the solutions 
obtained by these different solvers is practically indistin- 
guishable. Substantial savings in computer time can be 
realized by using the space-marching method in the 
PAB3D code if the analyses of many cases of jet flow 
conditions are required for design applications. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
May 23, 1996 
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Computational domain 


Inflow boundary 


Outflow boundary 


Turbulent shear mixing layer 
Jet centerline 



(a) Typical on-design jet flow configuration and terminology. 



(b) Single cell wedge grid for on-design jet flow computations. 

Figure 1. Sketch of typical axisymmetric on-design jet flow and computational grid. 



o 10 20 30 40 


x/R 

Figure 2. Centerline velocity decay for subsonic and supersonic jet flows computed with standard Jones-Launder k~E 
turbulence model with different compressibility corrections. 
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Figure 3. Axial velocity component distribution in cross sections in Mach 0.6 circular jet flow. 



Figure 4. Turbulence intensity distribution in cross sections in Mach 0.6 circular jet flow. 






Figure 5. Centerline velocity distribution for supersonic jet using k-£ turbulence model with different compressibility 
corrections. 



y/R 


Figure 6. Axial velocity component distribution along radial direction at x/R = 25 using k-£ turbulence model with 
different compressibility corrections. 
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Figure 7. Typical turbulence intensity distribution in axisymmetric jets computed by using Jones-Launder k-e turbu- 
lence model with Wilcox model of compressibility correction for u,/U e contours. 
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Mach number, M 

Figure 9. Potential core length as function of Mach number for subsonic and supersonic jets computed with Jones- 
Launder k-z turbulence model and different compressibility corrections. 
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Figure 10. Typical density contours and first shock-cell length definitions for circular underexpanded supersonic jet. Jet 
exit Mach number = 1 .50; pjp 0 - 1 .445. 
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Figure 1 1 . Computed centerline pressure distribution with different exit pressure ratios for Mach 2.0 jet computed with 
Jones-Launder k-E turbulence model. 



Figure 12. Computed and measured first shock-cell lengths for Mach 2.0 jet computed with Jones-Launder k~E turbu- 
lence model. 
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Figure 13. Centerline pressure distribution computed with Jones-Launder k-e model with different compressibility cor- 
rections for p e /p a = 1 .445. 



Figure 14. Centerline turbulence intensity computed with Jones-Launder k-e model with different compressibility cor- 
rections for p e !p a = 1 .445. 
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Figure 15. Space-marching and time-marching solutions for underexpanded Mach 2.0 supersonic jet computed with 
Jones-Launder k-e turbulence model. 
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(a) Mach number contours of typical underexpanded supersonic jet containing multiple Mach disks. 
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(b) Multiblock single cell wedge grid for jet flows containing multiple Mach disks. 

Figure 16. Flow configuration and computational grid for underexpanded jet containing one or more Mach disks. 
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(b) Mach number contour; Interval = 0.25. 

Figure 17. Density and Mach number contours for underexpanded circular jet containing multiple Mach disks 
Exit Mach number = 1.50; p e fp Q = 3.15; fine grid solution. 
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(a) Typical density contours of underexpanded elliptic jet in major plane of symmetry. 



(b) Adapted grid cross section at inflow plane. 



(c) Adapted grid longitudinal profile in plane of symmetry containing major axis of initial jet cross section. 
Figure 19. Adapted grid geometry for elliptic supersonic jet. Shape aspect ratio = 2.0. 
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(a) Plane of symmetry containing major axis; Interval = 0.20. 



(b) Plane of symmetry containing minor axis; Interval = 0.20. 

Figure 20. Mach number contours for underexpanded supersonic jet with elliptic exit cross section. Shape aspect 
ratio = 2.0; Exit Mach number = 2.00; NPR = 1 1.12; p e !p 0 = 1.445. 
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(a) Plane of symmetry containing major axis; Interval = 0.20. 



(b) Plane of symmetry containing minor axis; Interval = 0.20. 

Figure 21 . Mach number contours in major and minor planes of symmetry of elliptic jet shape ratio of 2.0; On-design 
NPR = 7.82. 
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(b) Centerline turbulence intensity distribution. 

Figure 28. Turbulence intensity distribution in circular jet computed by using single cell wedge grid. Exit Mach 
number = 2.0; on-design NPR = 7.82. 
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(b) Centerline turbulence intensity distribution. 

Figure 29. Turbulence intensity distribution in circular jet computed by using three-dimensional adaptive grids. Exit 
Mach number = 2.0; on-design NPR = 7.82 ; medium grid density: j, k = 40. 
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(b) Centerline turbulence intensity distribution. 

Figure 30. Turbulence intensity distribution in circular jet computed by using three-dimensional adaptive grids. Exit 
Mach number = 2.0; NPR = 7.82; high density adaptive grid: j 9 k = 56. 
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